Excitations in photoactive molecules from quantum Monte Carlo 



O 

o 

(N 
Sh 
Or 

< 

(N 



5-1 

43 



Friedemann Schautz 1 , Francesco Buda 2 , and Claudia Filippi 1 
1 Instituut-Lorentz, Universiteit Leiden, Niels Bohrweg 2, 2333 CA Leiden, The Netherlands 
2 Leiden Institute of Chemistry, Gorlaeus Laboratoria, 
P.O. Box 9502, 2300 RA Leiden, The Netherlands 
(Dated: February 2, 2008) 

Despite significant advances in electronic structure methods for the treatment of excited states, 
attaining an accurate description of the photoinduced processes in photoactive biomolecules is prov- 
ing very difficult. For the prototypical photosensitive molecules, formaldimine, formaldehyde and a 
minimal protonated Schiff base model of the retinal chromophore, we investigate the performance 
of various approaches generally considered promising for the computation of excited potential en- 
ergy surfaces. We show that quantum Monte Carlo can accurately estimate the excitation energies 
of the studied systems if one constructs carefully the trial wave function, including in most cases 
the reoptimization of its determinantal part within quantum Monte Carlo. While time-dependent 
density functional theory and quantum Monte Carlo are generally in reasonable agreement, they 
yield a qualitatively different description of the isomerization of the Schiff base model. Finally, we 
find that the restricted open shell Kohn-Sham method is at variance with quantum Monte Carlo 
in estimating the lowest-singlet excited state potential energy surface for low-symmetry molecular 
structures. 
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I. INTRODUCTION 

The absorption of visible light and its conversion to 
other forms of energy is at the heart of some of the most 
fundamental processes in biology. A familiar example 
of light absorption initiating a biological response is the 
primary event of vision: light induces a conformational 
change in rhodopsin, the photoreceptor in the retina, 
which is followed by a cascade of chemical reactions cul- 
minating in the stimulation of the optical nerve. A mi- 
croscopic understanding of light induced conformational 
changes in photoactive biomolecules is both important 
from a fundamental point of view and because of existing 
and potential applications in biology and biotechnology. 

The advances in understanding biological photosys- 
tems are so far mainly due to experimental discoveries 
since theoretical studies are currently hindered by the 
lack of a theoretical approach which is applicable to 
realistically large systems while possessing a sufficient 
degree of reliability. On the one hand, several accu- 
rate quantum chemical approaches have been developed 
for a proper description of excited states but they are 
only applicable to relatively small systems. For instance, 
complete active space second-order perturbation theory 
(CASPT2)i has been employed to investigate the pho- 
toisomcrization mechanism in simple models of the reti- 
nal chromophore of rhodopsin2ii4. The approach is able 
to accurately describe the excited state potential energy 
surface along the photoisomerization path, but it is lim- 
ited to relatively small model compounds and a proper 
description of the important ligand-protein interactions 
is still computationally prohibitive. On the other hand, 
density functional theory (DFT) based approaches have 
a much more favorable scaling with system size than 
CASPT2 and can therefore be applied to considerably 
larger molecules. In particular, the restricted open-shell 
Kohn-Sham method (ROKS) 5 ' 6 has been recently devel- 



oped to study the dynamics in low-spin excited states and 
used to model the full retinal chromophore, including rel- 
evant parts of the protein environment. The resulting 
excited state potential energy surface along the isomer- 
ization coordinate is qualitatively different from the one 
derived with the CASPT2 method^, though the model 
systems used in these two works are different and there- 
fore a direct comparison is not possible. Therefore, while 
the ROKS method is appealing for the low computational 
cost and for the possibility of performing molecular dy- 
namics in the excited state, its adequateness needs to be 
further validated. Alternatively, linear response calcu- 
lations within time-dependent density functional theory 
(TDDFT) 8 often yield accurate excitation energies but 
fail for instance in describing extended conjugated sys- 
tems^ or proton transfer— in excited states, that is, sys- 
tems closely related to photoactive molecules. The capa- 
bilities and limitations of TDDFT in describing excited 
state potential surfaces of conjugated organic molecules 
have been extensively investigated in Ref. 0. 

Quantum Monte Carlo (QMC) is an alternative to 
conventional quantum chemical and density functional 
methods, and has been successfully employed to compute 
ground state properties of large molecules and solidsi^. 
Compared to other theoretical approaches, QMC has the 
advantage that it can be applied to sufficiently large sys- 
tems and still provide an accurate description of both 
dynamical and static electronic correlation. Despite the 
successful use of QMC for ground state problems, there 
is relatively little experience on its application to excited 
utate^t^Hi^!^. The recent QMC computation of exci- 
tation energies of large silicon nanostructuresi 5 - is very 
encouraging but the simple HOMO-LUMO wave func- 
tions employed there are not likely to be adequate for 
photoactive systems due to the more complex nature of 
their electronic excitation. 

To compare the accuracy of ROKS, TDDFT and QMC 
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in the study of photochemical processes, we compute the 
excitation to the lowest singlet state for a set of proto- 
typical photoactive molecules: formaldimine (CH2NH), 
formaldehyde (CH2O) and a minimal protonated Schiff 
base model (CsHgNH^ ) of the retinal chromophore. For 
formaldimine and the protonated Schiff base model, wc 
find that ROKS differs quantitatively and qualitatively 
from the other methods under consideration at low- 
symmetry molecular structures. While TDDFT exci- 
tation energies are fairly accurate in most situations, 
this method gives a qualitatively different result along 
a complete-active-space self-consistent-field (CASSCF) 
minimum energy path for the isomerization of the pro- 
tonated Schiff base model. Finally, we find that QMC 
provides a reliable estimate of the lowest singlet exci- 
tation energies of the studied molecules, provided one 
makes an adequate and careful choice of the trial wave 
function. Although simple mean-field HOMO-LUMO 
Jastrow-Slater wave functions are not always adequate 
for these systems, we can recover accurate excitations 
energies by using a relatively small expansion in Slater 
determinants, whose orbitals and/or coefficients are re- 
optimized within QMC. 

In Sec. [HI we review the theoretical approaches em- 
ployed in this work. The computational details are 
given in Sec. IIIII and the numerical results are shown in 
Sec. llVAl and lTVBl Finally, in Sec. lTvTl we discuss the 
sensitivity of the QMC results to the choice of the trial 
wave function. 



II. THEORETICAL METHODS 

We briefly review the theoretical methods used in this 
work for the computation of excited states, and refer for 
more details to the literature. 

The restricted open-shell Kohn-Sham (ROKS) me- 
thod^ is a recent modification of the ASCF approach 
used for the computation of multiplet splittingsiiiSiiiLSii. 
In the ROKS approach, the energies of the states given by 
single determinants are not computed in separate calcula- 
tions as in ASCF, but the linear combination correspond- 
ing to the desired state of pure symmetry is directly min- 
imized under the constraint of orthogonality among the 
Kohn-Sham orbitals. In particular, the energy of an open 
shell singlet is estimated as E(s) = 2E(m) — E(t), where 
E(m) is the energy of the mixed singlet configuration, 
i.e. a single determinant having the open shell orbitals 
occupied with electrons of opposite spin, and E(t) the en- 
ergy of the corresponding triplet configuration. Within 
ROKS, the energy E(s) is optimized using conventional 
ground state density functionals and a common set of 
orthogonal orbitals is used for both contributions^. 

Both the ASCF and ROKS approaches offer a practi- 
cal recipe to the computation of excited states but they 
cannot be fully justified from a theoretical point of view 
and their validity must be empirically corroborated. An 
appealing feature of ROKS is that the method can be 



easily combined with ab-initio molecular dynamics and 
used to optimize the geometries in the excited state, ac- 
cess adiabatic excitations, and study the dynamics in the 
excited stato 5 i 22 i 23 i 24 . In general, even though the ROKS 
method tends to underestimate the excitation energies in 
particular for ir — > tt* transitions 22,25,26 , it was shown to 
give a good description of the optimal geometries of the 
lowest excited states of small organic molecules, espe- 
cially for n — > 7r* transition o 5 ! 22 . 

Time-dependent density functional theory (TDDFT) 
is a different framework for the calculations of excited 
state properties which has become widely used in re- 
cent years 8 -. The method can handle large systems 
and, differently from ASCF or ROKS, is formally ex- 
act even though, in practice, one has to resort to ap- 
proximate exchange-correlation functionals. TDDFT has 
been extensively applied to the computation of vertical 
excitation energies since the calculation of forces within 
TDDFT is not straightforward and only recently a few 
implementation and applications of TDDFT to compute 
excited state geometries and adiabatic excitations have 
been published 11 ! 27 ! 28 ' 29 . 

Several quantum chemical approaches have been devel- 
oped for a proper description of excited states. Methods 
such as multi-reference configuration interaction (MRCI) 
and complete active space second-order perturbation the- 
ory (CASPT2) rely on expanding, explicitly or implicitly, 
the wave function in Slater determinants. As the system 
size increases and the energies of the single-particle or- 
bitals become closely spaced, the space of orbitals which 
must be included in the expansion to recover a signifi- 
cant fraction of electronic correlation grows enormously. 
Therefore, these techniques are very accurate but can 
only be applied to small systems. Even though CASPT2 
was originally proposed as a method to compute excited 
state energies with an accuracy not better than 0.5 eV, it 
is now regarded as an approach which on average yields 
excitations in agreement with experiments to better than 
0.2 eVi The method is quite sensitive to the construc- 
tion of the active space which must include all important 
orbital excitations and is limited on current computers 
to a maximum of about 15 active orbitals. 

Quantum Monte Carlo techniques^ is an alternative 
to density functional and conventional quantum chem- 
istry approaches. While many studies have demonstrated 
the use and reliability of QMC for the description of 
ground state properties of molecular and solid systems, 
relatively little experience exists concerning its applica- 
tion to low-lying excited states. Recent studies of the ex- 
cited states of methane, ethene, and small hydrogenated 
Si clusters indicate that the method is capable of re- 
producing the excitation energies of accurate quantum 
chemistry calculationsiii^S. The QMC approach was 
also recently applied to the study of the excitations of 
large silicon nano-clusters, in combination with simple 
trial wave functions^ 5 -. QMC methods provide a stochas- 
tic solution of the Schrodinger equation: in diffusion 
Monte Carlo (DMC), the imaginary-time evolution op- 
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erator exp(— Ht) is used to project out the ground state 
from a given trial wave functional. To prevent the col- 
lapse to the bosonic ground state in fermionic systems, 
one works in the fixed-node approximation, that is, finds 
the best solution which has the same nodes as a given 
trial wave function. The solution is variational for the 
lowest state of a given spin symmetry belonging to a one- 
dimensional irreducible representation of the point group 
of the molecule. It is exact for any state if the nodes are 
exact. Therefore, if the nodal surface of the trial wave 
function is a good approximation to the excited state one, 
the fixed-node constraint can be used to access accurate 
excitation energies also of states which are not the lowest 
in their symmetry. 

The trial many-body wave function employed in this 
paper is of the Slater- Jastrow form: 

*T = ^2 d n D n D n II J ( r U > r « > »>) . 

n aij 

Djj and are Slater determinants of single particle or- 
bitals for the up- and down-spin electrons, respectively, 
and the orbitals are represented using atomic Gaussian 
basis. The Jastrow factor correlates pairs of electrons i 
and j with each other, and with every nucleus a, and dif- 
ferent Jastrow factors are used to describe the correlation 
with different types of atoms. The parameters in the Jas- 
trow factor are optimized within QMC using the variance 
minimization method 32 . The Jastrow factor is positive 
and does not alter the nodal surface of the wave function 
which is instead fixed by the determinantal parish 3 ". Par- 
ticular attention must therefore be paid to the choice of 
the Slater component which is usually a linear combina- 
tion of a small number of determinants. In the context of 
excited states, the complete-active-space self-consistcnt- 
field (CASSCF) variant of the multi-configuration self- 
consistent-field method (MCSCF) is particularly useful. 
These wave functions include all possible excitations for 
a given set of electrons within a chosen set of orbitals. 
When the excited state is not orthogonal to the ground 
state by symmetry, the determinantal component of the 
trial wave function is obtained in a state-average MC- 
SCF approach 3 ^, that is, by optimizing an average of the 
ground and excited state energies. Thus, the orbitals 
represent a compromise for describing both states. 

Since the optimal orbitals and expansion coefficients in 
the presence of the Jastrow factor may differ from their 
optimal values in its absence, it is important to reop- 
timize them in the presence of the Jastrow component. 
To this end, we extended the energy fluctuation poten- 
tial (EFP) method 35 to simultaneously minimize the en- 
ergy with respect to the orbitals and the expansion co- 
efficients of a Slater- Jastrow wave function, as well as 
to handle state averaging necessary for excited states^. 
In the absence of the Jastrow component, the method is 
analogous to the MCSCF technique for the lowest-state 
of a given symmetry, and to a state-average MCSCF ap- 
proach if the excited state of interest is not the lowest in 
its symmetry. Once the Jastrow factor is included, the 



orthogonality between the ground and excited states is 
only approximately preserved in the state-average EFP 
approach. The approach was tested for several singlet 
states of ethene and was shown to systematically improve 
the starting trial wave functions, correcting the initial ex- 
citation energies by as much as 0.5-0.6 eV and yielding 
results in excellent agreement with experiments^. 



III. COMPUTATIONAL DETAILS 

The ground-state DFT, and the excited-state ROKS 
and TDDFT calculations are performed with the Car- 
Parrinello molecular dynamics CPMD cod o 36 i 37 . Wc 
employ the BLYP generalized gradient approximation 
for the exchange and correlation functionate* 3 ^, the 
Goedecker pseudopotentials-^, an energy cutoff of 70 Ry 
for the plane-wave expansion, and a box size about 5 A 
larger then the size of the molecule. In order to avoid the 
inherent periodicity of a plane-wave calculation, we use 
the method described in Ref4i, which solves the Pois- 
son equation for non-periodic boundary conditions, thus 
enabling the study of isolated molecules. 

For formaldimine, the multi-reference configuration 
interaction singles and doubles (MR-CISD) calcula- 
tions and the optimization of the excited state ge- 
ometry within the state-average CASSCF method are 
performed with the COLUMBUS quantum chemistry pro- 
gram^. Equal weights are used in the state-average 
CASSCF calculations for the optimization of the geome- 
tries. The reference space for MRCI is of 6 active elec- 
trons in 6 orbitals and the final MRCI energetics in- 
clude Davidson corrections. It must be stressed that 
these MRCI calculations were performed with a moderate 
basis ((10 s6p3d)/[4s3pld] for carbon and nitrogen, and 
(7s3p)/[2slp\ for hydrogen) and could certainly be im- 
proved. However, for the purpose of establishing the re- 
liability of the other theoretical approaches, we consider 
the accuracy of the MRCI energetics to be sufficient. 

For the QMC calculations, we use the CHAMP quan- 
tum Monte Carlo code^ 3 - and norm-conserving sp-non- 
local pseudopotentials for carbon, nitrogen and oxigen, 
generated in an all-electron Hartree-Fock calculation for 
the atoms^i. The orbitals in the determinantal compo- 
nent of the wave functions are expanded in the Gaus- 
sian basis sets (llsllp2d) /[4s4p2d] for carbon, nitrogen, 
and oxigen, and (10s2p)/[3s2p\ for hydrogen. The basis 
sets are optimized at the HF level for formaldimine and 
formaldehyde. The determinantal part of the wave func- 
tion, before reoptimization in QMC, is generated within 
Hartree-Fock, CASSCF or state-average CASSCF, using 
the quantum chemistry package GAMESS(US) 45 . Equal 
weights are used in the state-average CASSCF calcula- 
tions, and in the state-average EFP optimization of the 
wave function. The Jastrow factor contains electron- 
electron, electron-nucleus and electron-electron-nucleus 
terms and is described in Ref. ^ij For reasons of ef- 
ficiency, most calculations are performed omitting the 
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electron-electron-nucleus terms since the excitation en- 
ergies for these systems computed with or without the 
three-body terms are the same within better than 0.1 
eV—. The diffusion Monte Carlo time-step used for these 
molecules is 0.075 H _1 . Most of the QMC results pre- 
sented below are obtained in diffusion Monte Carlo. Vari- 
ational Monte Carlo (VMC) is also used to compute var- 
ious expectation values of the trial Jastrow-Slater wave 
function. 



IV. RESULTS 

The photosensitive molecules we investigate are 
schematically shown in Fig. ^ In formaldimine and 
formaldehyde, the lowest singlet excitation has predomi- 
nantly a n — > 7r* character and, in the protonated Schiff 
base model, a n — > ir* character. The performance of the 
DFT-based approaches may differ for the two types of 
excitation, as has previously been stated for the ROKS 
method. 

While QMC does not seem to be sensitive to the char- 
acter of the excitation, a different complication is encoun- 
tered when performing excited state QMC calculations. 
If the excited state of interest is the lowest state of a 
given spin symmetry belonging to a one-dimensional ir- 
reducible representation, the DMC energy is variational. 
In all other cases, DMC is no longer variational and 
the quality of the trial wave function becomes increas- 
ingly important. The vertical and adiabatic excitations 
of formaldimine and formaldehyde belong to the first cat- 
egory while the excitations of the minimal protonated 
Schiff base model and of formaldimine along its isomer- 
ization path belong to the other case. 
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FIG. 1: Structure of the investigated molecules. In 
formaldimine and the protonated Schiff base model, the iso- 
merization is around the bond indicated with an arrow. 
H1CNH is the dyhedral angle varied in formaldimine. 



A. Formaldimine and formaldehyde 

In the n — > ir* excitation of formaldimine and 
formaldehyde, a lone-pair electron is transferred to a 



7r* antibonding orbital. The excitation is almost purely 
of the HOMO-LUMO type and has therefore been con- 
sidered ideal for the ROKS approach-, which was also 
used to study the excited state cis-trans isomerization of 
formaldimine in a Born-Oppcnheimcr molecular dynam- 
ics simulation 5 - and, more recently, in a non-adiabatic 
Car-Parrinello dynamics^. 



TABLE I: Vertical and adiabatic lowest singlet excitation en- 
ergies in eV for formaldehyde and formaldimine, calculated 
within ROKS, TDDFT and DMC. The numbers in parenthe- 
ses are the statistical errors on the DMC results. 



system excitation ROKS TDDFT 


DMC 


Expt 


CH 2 vertical 3.58 


3.90 


4.24(2) 3.94' 


\ 4.07 i \ 4.2 C 


adiabatic 3.13 


3.51 


3.74(2) 


3.50 d 


CH 2 NH vertical 4.63 


5.34 


5.32(2) 


5.0-5.4 6 


adiabatic 2.85 


3.23 


3.21(2) 





Ref . El b Ref . |H c Ref . |U d Ref. HI 

In TableHJ we list the vertical and adiabatic lowest sin- 
glet excitation energies, evaluated using ROKS, TDDFT 
and DMC. The vertical excitations are computed on the 
ground state DFT geometries, while the adiabatic exci- 
tations on the geometries optimized in the excited state 
using ROKS. The adiabatic geometry of formaldehyde is 
known experimentally and is well reproduced by ROKS 5 . 
Vertical and adiabatic transitions are underestimated by 
ROKS by as much as 0.5 eV, while the TDDFT results 
are in reasonable agreement with experiments. These 
findings are consistent with previous ROKS calculations 
for both molecules^, and with TDDFT calculations of the 
vertical 5 "^ and adiabatic^ excitations of formaldehyde. 

The DMC excitations are obtained using a compara- 
ble description of the ground and excited states. A one- 
determinant trial wave function is used for the ground 
state, and a two-determinant singlet wave function for 
the excited state, corresponding to a single excitation 
from the doubly-occupied n HOMO to the ir* LUMO. 
The starting orbitals in the determinantal component of 
the QMC wave function are from a HF calculation in 
the ground state, and a two-determinant MCSCF cal- 
culation in the excited state. For both states, all or- 
bitals are subsequently optimized in the presence of the 
Jastrow factor with the EFP method. For formalde- 
hyde, the DMC excitation energies are slightly higher 
than available experimental numbers and results from 
highly-correlated quantum chemistry calculations, which 
however show a significant spread. The vertical excita- 
tion energies computed with quantum chemistry tech- 
niquea 52 i 53 i 54 i 55 range between 3.98 eV from EOM-CC 
and 4.19 eV from MRCI^ 5 ., while MRCI calculations for 
the adiabatic transition 5 ^ yield an excitation energy of 
3.60-3.66 eV. For formaldimine, the DMC vertical and 
adiabatic excitations are in good agreement with MRCI 
calculations 5 ^. 

While the success of DMC in describing these vertical 
and adiabatic excitations is encouraging, it is important 
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to assess its performance when variationality is lost as 
happens along the low-symmetry isomerization path in- 
duced by the excitation. We therefore consider the proto- 
typical case of the isomerization of formaldimine around 
the C-N double bond. The isomerization path is con- 
structed by constraining the torsional angle H1CNH (see 
Fig. at values between and 90 degrees, with incre- 
ments of 15 degrees. The molecule has C s symmetry at 
and 90 degrees, and no symmetry at intermediate angles. 



Formaldimine: ROKS geometries 




'0 15 30 45 60 75 90 

Torsional angle (deg) 

FIG. 2: Lowest-singlet excitation energies of formaldimine 
in eV calculated with ROKS, TDDFT, MRCI and DMC on 
the excited state geometries optimized with ROKS at con- 
strained torsional angles. The excitations computed within 
a two-determinant MCSCF calculation (jj. — ||) are also 
shown. The statistical error on the DMC results is smaller 
than the size of the symbols. 

In Fig. El we show the ROKS, TDDFT, DMC and 
MRCI excitation energies on the excited state geome- 
tries optimized with ROKS at constrained torsional an- 
gles. The excitation energies are given with respect to 
the ground state energy consistently computed within 
the same approach on the DFT ground state geometry at 
zero torsional angle. The DMC excited state energies are 
obtained with a trial wave function from a state-average 
CASSCF with an active space of 6 electrons in 6 orbitals, 
whose expansion coefficients are then reoptimized in the 
presence of the Jastrow factor with a state-average EFP 
method. The DMC ground state energy at zero torsional 
angle is computed with an unoptimizcd HF dcterminan- 
tal component. The DMC excitations are in very good 
agreement with the MRCI values, with a maximum de- 
viation of 0.13 eV along the curve. 

While the TDDFT excitations agree with the MRCI 
values to better than 0.2 eV, the ROKS curve differs sig- 
nificantly. In particular, MRCI gives a barrier to isomer- 
ization along the geometries corresponding to an energy 
minimum path in ROKS. One can possibly understand 
the behavior of ROKS by looking at the results obtained 
with a two-determinant MCSCF (without state-average) 
approach along the same path. As shown in Fig. El the 
two-determinant MCSCF curve is qualitatively very sim- 



ilar to the ROKS curve. For the two-determinant MC- 
SCF calculation, only the orthogonality constraint on the 
open shell orbitals keeps the wave function from com- 
pletely collapsing to the ground state. By analogy, the 
ROKS approach is likely to suffer from the same prob- 
lem whenever ground and excited states do not belong to 
different irreducible representations 22 . 



Formaldimine: CASSCF geometries 
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FIG. 3: Lowest-singlet excitation energies of formaldimine 
in eV calculated with ROKS, TDDFT, MRCI and DMC on 
the excited state geometries optimized using a state-average 
CASSCF (see text) at constrained torsional angles. 

To further investigate the constraint isomerization 
path of formaldimine, we optimize the geometries using 
the excited-state forces from a state-average CASSCF 
approach with an active space of 6 electrons in 6 or- 
bitals. As already pointed out in early MRCI studies 
by Bonacic-Koutecky et al£^, to properly describe the 
isomerization of formaldimine, one should map the po- 
tential energy surface with respect to the CNH valence 
angle and a properly symmetrized dyhedral angle. How- 
ever, the path obtained within CASSCF by only con- 
straining the H1CNH dyhedral angle is reasonably close 
to the optimal path. We find that the main difference be- 
tween the ROKS and CASSCF paths is in the behavior 
of the angle CNH which, in ROKS, takes his final value 
corresponding to a torsional angle of 90 degrees as soon 
as the molecule is displaced from planarity. 

The excitations computed with TDDFT, ROKS, DMC 
and MRCI on the CASSCF geometries are shown in 
Fig. El The DMC calculations are performed with the 
same type of wave function previously used for the ROKS 
path. The energy barrier to isomerization present in 
Fig. disappears in MRCI as this barrier was an ar- 
tifact of using the geometries optimized within ROKS. 
The DMC excitation energies are very close to the MRCI 
values with a maximum difference of 0.1 eV along the 
CASSCF path. TDDFT is in reasonable agreement with 
QMC also along this path. For the CASSCF geometries, 
ROKS calculations produce a curve of similar shape as 
those obtained with the other methods, but significantly 
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shifted toward lower energies. 



B. Protonated Shiff base model 

The C5H6NHJ protonated Shiff base molecule repre- 
sents a minimal model for studying the retinal photoiso- 
merization process in rhodopsin. Given its relevance and 
combined simplicity, this molecule is ideal for accessing 
the relative accuracy of different theoretical approaches. 
Moreover, this model has been extensively studied within 
CASPT2 using geometries optimized in the excited state 
with CASSCBM and, more recently, with CASPT21 

Since ROKS was previously employed to study the ex- 
cited state of the full retinal chromophorc including rele- 
vant parts of the protein environment, it is interesting to 
use the same approach to optimize the structure of this 
simpler model. In Fig. g| we show the ROKS, TDDFT 
and DMC energetics computed on the geometries opti- 
mized within ROKS along the relevant isomerization co- 
ordinate represented by the torsional angle around the 
central C-C double bond (see Fig. When optimiz- 
ing the excited state geometry with ROKS, the molecule 
remains planar and the main effect of the excitation is 
a considerable lengthening of the double bonds and a 
shortening of the single bonds, thus reversing the con- 
jugation of the molecule. The ROKS potential energy 
surface along the torsion is quite flat with a maximum 
at 90 degrees. This behavior is qualitatively different 
from the CASSCF and CASPT2 energy profile^, where 
the torsion accelerates the system towards the conical 
intersection, thus spontaneously inducing the photoiso- 
merization. Therefore, while the ROKS method shows a 
stretching mode starting from the Franck-Condon region 
similar to the CASSCF result, it does not reproduce the 
qualitative shape of the excited state CASSCF potential 
energy surface along the torsional mode. 

The DMC excited state energies in Fig. 0] are com- 
puted on the ROKS geometries with a trial wave function 
from a state-average CASSCF with an active space of 6 
electrons in 6 orbitals, whose expansion coefficients are 
then reoptimized in the presence of the Jastrow factor 
with a state-average EFP method. The TDDFT exci- 
tation energies are higher than the ROKS values by as 
much as 2 eV, and in agreement with the DMC results 
to better than 0.2 eV. The TDDFT and DMC poten- 
tial energy curves have a very different shape than the 
one obtained within ROKS. In the protonated Schiff base 
model, the ground and excited states belong to the same 
irreducible representation both when the molecule is pla- 
nar and twisted. The behavior of ROKS can possibly 
be explained as due to a contamination of the excited 
state with the ground state as in the case of twisted 
formaldimine. 

To allow for a comparison with existing CASPT2 cal- 
culations on this model, we consider three geometries 
which were optimized in Ref. |2j within state-average 
CASSCF and where the CASPT2 energies are also avail- 
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FIG. 4: Lowest-singlet excitations energies for the protonated 
Schiff base model in eV calculated with ROKS, TDDFT and 
DMC on the excited state geometries optimized with ROKS 
at constrained torsional angles. The excitation energies are 
given with respect to the ground state energy consistently 
computed within the same approach on the DFT ground state 
cis-geometry at zero torsional angle. 

able. These structures correspond to the ground state cis- 
configuration where the Franck-Condon (FC) excitation 
is computed, to the geometry which demarcates where 
torsion becomes dominant along the isomerization path 
(denoted with HM in Ref.0), and to the So/Si conical in- 
tersection (CI). Without a direct comparison with exper- 
imental data, it is difficult to access the accuracy of these 
excited state structures: for instance, when compared to 
geometries optmized with CASPT2, the CASSCF struc- 
tures are very similar at the conical intersection but sig- 
nificantly different at constrained planar symmetry^. 

TABLE II: Lowest-singlet excitation energies for the pro- 
tonated Schiff base model in eV, calculated with TDDFT, 
CASPT2 and DMC on the ground state cis-configuration 
(FC), on the geometry (HM) which demarcates where tor- 
sion becomes dominant, and on the conical intersection (CI). 
The CASSCF geometries and the CASPT2 numbers are from 
Ref. y. The excitation energies are given with respect to the 
ground state energy consistently computed within the same 
approach on the CASSCF ground state cis-geometry at zero 
torsional angle. 

Geometry TDDFT CASPT2 DMC 
~FC 3?90 4~02 4.38(5) 

HM 4.12 3.71 4.22(5) 

CI 2.18 2.19 2.58(5) 

In Table CH we list the TDDFT, CASPT2 and DMC 
excitation energies at the FC, HM and CI geometries. 
The DMC calculations are performed with the same type 
of wave function previously used for the ROKS path. The 
use of larger active spaces (6 electrons in 9 orbitals or 8 
electrons in 8 orbitals) and the reoptimization of the ac- 
tive orbitals with the state-average EFP method yield 
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DMC energies compatible to better than 0.1 eV. While 
the CASPT2 and QMC results are qualitatively similar, 
the CASPT2 energies are lower than the QMC values by 
as much as 0.5 eV. The order of the TDDFT excitation 
energies at the FC and HM configurations are instead 
reversed with respect to the DMC values: the TDDFT 
excitation is lower at FC than at HM, so TDDFT gives 
a barrier to isomerization along the CASSCF path. A 
valid question is whether this barrier survives when us- 
ing an excited state path fully optimized within TDDFT. 
Recently, it has been shown that the TDDFT gradient 
for various protonated Schiff base models differs qualita- 
tively from that of CASSCF/CASPT2, driving the sys- 
tem from the FC point to a planar fictitious stationary 
ponrfcii. 



Protonated Schiff base model: twist geometries 

4.5 - 




Torsional angle (deg) 

FIG. 5: Excitation energies for the protonated Schiff base 
model in eV, calculated with TDDFT and DMC on a set of 
geometries generated by rigidly increasing the torsional an- 
gle, from the HM configuration. The TDDFT, DMC and 
CASPT22 energies at FC and HM are also given. 

Finally, in order to further compare TDDFT and 
QMC, we generate a set of geometries for CsHeNH^" by 
starting from the HM structure of Ref. and increasing 
the torsional angle up to about 90 degrees while keep- 
ing all the other internal coordinates fixed. In Fig. [SJ we 
show the TDDFT and DMC energies, and the CASPT2 
results at FC and HM. Along the torsional path after 
HM, TDDFT and DMC follow closely each with a larger 
deviation at the end of the path. 



C. Sensitivity of DMC to the trial wave function 

Using as examples the vertical excited state and the 
adiabatic isomerization path of formaldimine, we demon- 
strate how sensitive the QMC energies are to the choice 
of the wave function and how this sensitivity can vary 
along the excited state potential energy curve. 

The vertical lowest-singlet excited state of formaldimi- 
ne does not have a strong multi-configurational charac- 



ter, and a two-determinant Jastrow-Slater wave function 
to preserve spin symmetry is found to be sufficient for 
this particular state. The QMC energies are variational 
since this excited state is the lowest in its symmetry, and 
orthogonality between ground and excited state is au- 
tomatically ensured. For the ground state, a single de- 
terminant wave function gives an adequate description. 
In Table GUI we show the VMC and fixed-node DMC 
energies determined with different choices of orbitals in 
the determinantal component of the wave function. The 
starting trial wave function uses orbitals obtained from a 
HF and a two-determinant MCSCF calculations for the 
ground and excited state, respectively. By optimizing the 
orbitals with the EFP method, the VMC energy drops 
by 10 mHartree in the ground state and by 15 mHartree 
in the excited state. However, the gain in the DMC en- 
ergies is only of a few mHartree and is actually more 
significant in the ground state. The resulting DMC ex- 
citation energy is only slightly higher as a result of the 
optimization. 

TABLE III: VMC and DMC ground state (So) and lowest- 
singlet excited state (Si) energies in Hartree for formaldimine, 
calculated at the ground state geometry. In the Jastrow- 
Slater wave function, a single determinant is used for the 
ground state and two determinants for the excited state. The 
DMC excitation energies in eV are computed using unopti- 
mized (HF for S and MCSCF for Si) and optimized (EFP) 
orbitals for both states. 



State Orbitals 


-Evmc 


-Bdmc 


AE (eV) 


So HF 


-17.2973(4) 


-17.3685(5) 




EFP 


-17.3082(4) 


-17.3726(5) 




Si MCSCF 


-17.1185(4) 


-17.1756(5) 


5.25(2) 


EFP 


-17.1334(4) 


-17.1772(4) 


5.32(2) 



Along the isomerization path of formaldimine, orthog- 
onality between ground and excited state is no longer 
maintained and a higher sensitivity of the QMC results 
to the trial wave function may be expected than in the 
case of the vertical excitation. In Fig. we compare 
the DMC excitation energies along the ROKS isomeriza- 
tion path of formaldimine for different choices of wave 
functions previously employed in other QMC studies of 
excited states. At and 90° torsional angles where the 
energy is variational due to symmetry, the spread of the 
DMC energies due to the use of different wave functions is 
significantly smaller than at intermediate angles. A sim- 
ple two-determinant HOMO-LUMO wave function with 
HF orbitals shows a discrepancy as large as 1.5 eV with 
our best DMC results obtained with a 6 electrons in 6 or- 
bitals CASSCF wave function whose CI coefficients have 
been reoptimized with the state-average EFP method. 
The wave function denoted with CIS1 includes all sin- 
gle excitations from the HOMO, and can be resummed 
to two determinants, where only the LUMO has there- 
fore been changed with respect to the HF orbitals. The 
CIS1 energies represent an improvement at the end points 
of the path but remain as poor as those obtained with 
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a HOMO-LUMO wave function at almost all other an- 
gles. If all single excitations are included in a CIS wave 
function, the excitation energies are significantly closer to 
the CASSCF-EFP results along the whole path, with an 
almost constant discrepancy of 0.3-0.5 eV. Finally, one 
could be tempted to use a two-determinant wave func- 
tion obtained in a MCSCF calculation (without state- 
average). While this wave function performs well at 
and 90 degrees where ground and excited states are or- 
thogonal by symmetry, it represents a poor starting point 
at low-symmetry configurations as already discussed in 
Section llV Al yielding DMC energies which are obviously 
non variational. 



Formaldimine: DMC energies 

Wave function: * HF HOMO-LUMO 

v CIS1 
► CIS 

* MCSCF U-iT 

• CASSCF-EFP 
□ MRCI 




30 45 60 

Torsional angle (deg) 



FIG. 6: DMC lowest-singlet excited state energies of 
formaldimine in eV, computed on the ROKS geometries at 
various torsional angles, using different trial wave functions. 
See text for more details. 

Finally, the effect of truncating the determinantal ex- 
pansion according to a threshold on the coefficients is 
investigated. It is indeed customary in QMC to apply 
a threshold for computational efficiency, justified by the 
very different role of the reference wave function in QMC 
compared to conventional quantum chemistry methods. 
A smaller number of determinants is needed in a Jastrow- 
Slater wave function since the reference wave function 
does not define the single-particle excitation space for 
the description of dynamical correlation as is the case for 
a method like MRCI. Moreover, one hopes that the ef- 
fect of determinants with a small coefficient on the nodal 
surface of the total wave function is not significant. 

In Table |TV| we show the VMC and DMC excited 
state energies for formaldimine, computed on the ROKS 
geometries at various torsional angles when applying 
two different thresholds on the expansion coefficients in 
symmetry-adapted configuration state functions. The 
starting trial wave function is obtained from a state- 
average CASSCF with an active space of 6 electrons in 
6 orbitals. As the threshold is lowered from 0.1 to 0.01, 
both VMC and DMC energies become higher at all an- 
gles. Since at and 90 degrees the energies are varia- 
tional due to symmetry, one is unequivocally aiming at 



TABLE IV: VMC and DMC lowest-singlet excited state ener- 
gies for formaldimine, computed on the ROKS geometries at 
various torsional angles. Different determinantal components 
are used in the trial wave functions, with thresholds of 0.1 
and 0.01 on the expansion in symmetry-adapted configura- 
tion state functions from a state-average CASSCF, and with 
CASSCF and EFP-optimized expansion coefficients. 



Threshold 


0.1 


0.01 


0.01 


Coefficients 


CASSCF 


CASSCF 


EFP 


Angle (deg) 


Number of determinants 





4 


42 


23 


30 


9 


132 


46 


60 


8 


108 


54 


90 


4 


71 


35 


Angle (deg) 


VMC 


energies (Hartree) 





-17.158(1) 


-17.152(1) 


-17.165(1) 


30 


-17.149(1) 


-17.144(1) 


-17.158(1) 


60 


-17.180(1) 


-17.178(1) 


-17.190(1) 


90 


-17.200(1) 


-17.193(1) 


-17.205(1) 


Angle (deg) 


DMC energies (Hartree) 





-17.2099(5) 


-17.2063(5) 


-17.2113(4) 


30 


-17.2027(5) 


-17.2000(5) 


-17.2062(4) 


60 


-17.2338(5) 


-17.2313(5) 


-17.2360(4) 


90 


-17.2502(4) 


-17.2474(5) 


-17.2527(4) 



obtaining the lowest possible energy at those geometries 
and one would have expected a lowering of the energy by 
including more configurations. This indicates that the 
result is strongly dependent on the chosen threshold if 
one does not reoptimize the determinantal expansion in 
the presence of the Jastrow factor. The coefficients of 
the starting CASSCF wave function are therefore reop- 
timized with the state-average EFP method. The natu- 
ral orbitals of the averaged single-particle density matrix 
of the reoptimized expansions are here used to obtain a 
more compact wave function, and a threshold of 0.01 is 
then applied. The corresponding VMC and DMC ener- 
gies are also shown in Table HVl At all angles, the VMC 
energies for the reoptimized wave function are lower than 
the values obtained using the original CASSCF coeffi- 
cients with respect to the same threshold. Moreover, the 
optimal energies are also systematically better than the 
VMC values obtained with a threshold of 0.1. In Ta- 
ble ]W\ we also list the number of determinants with co- 
efficients greater than the chosen threshold. As expected, 
due to the inclusion of dynamical correlation through the 
Jastrow factor, the wave function becomes more compact 
as an effect of the reoptimization. The DMC energies be- 
have similarly to the VMC values with respect to both 
threshold and reoptimization. The excitation energies 
obtained in DMC with the reoptimized wave function are 
in excellent agreement with the MRCI values as shown 
in Section TlV Al If a threshold of 0.1 is used when re- 
optimizing the expansion coefficients in a state-average 
EFP method, there is no improvement in the QMC en- 
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ergies compared to the values obtained with the original 
CASSCF coefficients and the same threshold. 

Finally, if the orbitals are optimized with the state- 
average EFP approach and a threshold of 0.1, both VMC 
and DMC energies improve and became equal to the val- 
ues obtained with the CASSCF-EFP with 0.01 threshold. 
For instance, for a torsional angle of 30 degrees, the opti- 
mization of the orbitals yields a VMC and a DMC energy 
of -17.156(1) and -17.2071(4) Hartree, respectively. We 
want to stress that there is in general no justification for 
using a threshold as high as 0.1 and the apparent agree- 
ment with the optimized energies is here a fortunate case. 

V. CONCLUSIONS 

Using TDDFT, ROKS and QMC, we have investigated 
the lowest-singlet excitation energies along various iso- 
merization paths for the following representative pho- 
toactive molecules: formaldehyde, formaldimine and a 
minimal protonated Schiff base model CsHgNH^. 

We show that fixed-node diffusion Monte Carlo can 
give accurate excitation energies, provided a careful 
choice of QMC trial wave function is made. While simple 
HOMO-LUMO trial wave functions are not always ade- 
quate to describe the excited states of these photoactive 
molecules, accurate results are recovered when using a 
relatively small expansion in Slater determinants, whose 
coefficients and/or orbitals are reoptimized in the pres- 
ence of the Jastrow factor with the EFP method. 



TDDFT yields excitation energies which are generally 
in reasonable agreement with the QMC results. However, 
the TDDFT energies for the minimal model of the retinal 
chromophore are in qualitative disagreement with QMC 
and CASPT2, giving a barrier to isomerization along the 
CASSCF minimal energy path. 

We find that the ROKS method does not produce reli- 
able results for the excited-state potential energy surface 
at low-symmetry configurations. The major source of er- 
ror in the ROKS approach seems to be the contamination 
of the excited state with the ground state. For example, 
ROKS predicts an energy barrier to isomerization with 
a maximum at 90 degrees along the relevant torsional 
angle of the minimal protonated Schiff base model of 
the retinal chromophore, while TDDFT and QMC show 
a minimum at this point. Therefore, even though the 
ROKS method is appealing for its simplicity in comput- 
ing forces, it should be generally used with caution in 
excited-state molecular dynamics simulations. 
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